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ABSTRACT 

We use a set of hydrodynamical (Hydro) and dark matter only (DMonly) simulations 
to calibrate the halo mass function (HMF). We explore the impact of baryons, propose an 
improved parametrization for spherical overdensity masses and identify differences between 
our DMonly HMF and previously published HMFs. We use the Magneticum simulations, 
which are well suited because of their accurate treatment of baryons, high resolution, and 
large cosmological volumes of up to (3818 Mpc) 3 . Baryonic effects globally decrease the 
masses of galaxy clusters, which, at a given mass, results in a decrease of their number density. 
This effect vanishes at high redshift z ~ 2 and for high masses M 2 oom > 1O 14 M0. We 
perform cosmological analyses of three idealized approximations to the cluster surveys by 
the South Pole Telescope (SPT), Planck , and eROSITA. We pursue two main questions: (1) 
What is the impact of baryons? - For the SPT-like and the Planck -like samples, the impact of 
baryons on cosmological results is negligible. In the eROSITA-like case, however, neglecting 
the baryonic impact leads to an underestimate of fl m by about 0.01, which is comparable 
to the expected uncertainty from eROSITA. (2) How does our DMonly HMF compare with 
previous work? - For the Planck -like sample, results obtained using our DMonly HMF are 
shifted by A (as) — A(as(fl m /0.27) 0 3 ) ~ 0.02 with respect to results obtained using the 
Tinker et al. (2008) fit. This suggests that using our HMF would shift results from Planck 
clusters toward better agreement with CMB anisotropy measurements. Finally, we discuss 
biases that can be introduced through inadequate HMF parametrizations that introduce false 
cosmological sensitivity. 

Key words: methods: numerical - galaxies: clusters: general - cosmological parameters - 
cosmology: theory - dark matter - large-scale structure of the Universe 


1 INTRODUCTION 

Galaxy clusters are the largest collapsed objects in the Universe. 
Their distribution in mass and redshift is highly sensitive to key 
cosmological parameters such as the matter density S2 m , or the 
amount of matter fluctuations in the Universe as (e.g. (Henry &| 
|Amaud|1991[|White, Efstathiou & Frenk|1993[ >. Furthermore, they 
can be used to constrain models of dark energy, the cosmic growth 


rate, and the neutrino sector fWang & Steinhardt]|1998[ Haiman, 
|Mohr & Holder|2001|. Catalogues from different cluster surveys 
have proven to be useful cosmological probes (e.g.|Vikhlinin et al.] 


2009, Mantz et al. 2010; Rozo et al.||2010; Benson et al. 2013 

Hasselfield et al. 2013 

Planck Collaboration et al. 2015b, Bocquet 

et al. 2015; Mantz et a 

L 20151. 


* E-mail: bocquet@usm.lmu.de 


The predicted abundance of galaxy clusters is linked to the lin¬ 
ear matter power spectrum through the halo mass function (HMF), 
which was first estimated analytically i jPress & Schechter||l974| l. 
Since then, numerical TV-body simulations have been used to cal¬ 
ibrate fitting functions (e.g. Jenkins et al. 2001 |Sheth, Mo~& Tor- 


]|2001^ White, Hernquist & Springel 2002 JReedeMdJ|2003 


Warren et a 


2008 

Crocce et al. 2010 

2011 

Angulo et al. 2012 


studies focus on the friends-of-friends (FoF) halo definition (Davis 
|et al.|1985) . However, real cluster samples are typically defined in 
terms of spherical overdensity (SO) masses. Only very few HMFs 
exist for different overdensity definitions (e.g. |Tinker et al.|2008~j 
|Watson et al.|2013| hereafter Tinker08, Watsonl3), and the former 
has developed into the standard reference used in most cluster cos¬ 
mology analyses. 
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For a particular FIMF parametrization to be useful in cosmo¬ 
logical studies, it is crucially important that it correctly captures 
both the variation in the HMF with redshift, and the sensitivity to 
cosmological parameters of interest such as the matter density fl m , 
the dark energy density IJa, the dark energy equation of state pa¬ 
rameter w, and ag. An ideal situation would be for the FIMF shape 
parametrization to be universal, where the variation with cosmol¬ 
ogy would be entirely captured by the cosmological sensitivity of 
the linear power spectrum of density fluctuations. For a FoF halo 
definition with linking length b = 0.2, or SO Aisomean, the FIMF 
has been found to be approximately universal over a wide range 
of redshifts and cosmologies ( [Jenkins et al.|200T ). More recently, 
|Bhattacharya et al.j {2011) have used a set of ACDM and tuCDM 
simulations to investigate the dependence of the FoF FIMF on cos¬ 
mological parameters. Their fit is accurate to 2% for ACDM, and 
it describes the wCDM HMF to within 10%. Similar results have 
also been reported in |Courtin et al.||201l| , although with slightly 
larger uncertainties. 

Any HMF obtained from iV-body dark matter only simula¬ 
tions potentially suffers from some bias because the baryonic com¬ 
ponent of the clusters is neglected. Recently, various authors have 
investigated the baryonic impact on the halo HMF using hydro¬ 


dynamic simulations (e.g. 

Cui et al. 

2012; Cui, Borgani & Mu- 

rante|2014 

Cusworth et a 

. 2014; M 

artizzi et al. 2014; Velliscig 

et al. 2014 

Vogelsberger et al. 2014 

Schaller et al. 2015). Their 


conclusions are highly sensitive to the details of the treatment of 
the baryonic component. For example, models without feedback 
from active galactic nuclei (AGN) lead to higher cluster masses 
(or higher abundance at fixed mass) than dark matter only simu¬ 
lations. Adding AGN feedback, however, leads to HMFs that are 
up to 20% lower than, or about equal to the dark matter only pre¬ 
diction. Also, these baryonic effects are stronger for low cluster 
masses, and hardly affect the high-mass population. These shifts in 
the predicted HMFs nearly correspond to the level of uncertainty 
from current cluster abundance measurements. Therefore, studies 
of the baryonic impact on the HMF are extremely important for 
progress in cluster cosmology. 

In this work, we use haloes extracted from the hydrody- 
namical Magneticum simulations (Dolag et al., in prep.; see also 
Hirschmann et al.|2014[|Saro et al.|2014[|Dolag, Komatsu & Sun-| 


yaev|20 15 Teklu et al. 2015). In combination with their dark matter 
only counterparts, these simulations allow us to investigate several 
key aspects of the HMF. In particular, we focus on the impact of 
baryons, universality for various SO definitions, and a compari¬ 
son of our dark matter only HMF with previously published fits. 
In Section [2] we present the Magneticum simulations and describe 
how the cluster samples are extracted. We discuss theoretical as¬ 
pects of the HMF in Section [3] where we also present a novel ap¬ 
proach for parametrizing the HMF for SO different from A 2 oom- 
In Section [4] we present our HMFs for different SO definitions, 
discuss the baryon impact and investigate differences between our 
dark matter only HMF and previous work. The cosmological im¬ 
pact of all these effects is presented in Section[5] We summarize our 
results in Section[6] where we also present step-by-step instructions 
on how to use our HMF fitting formulae. 

We consider the following SO definitions: (1) mean overden¬ 
sity mass M 200111 , which is the mass enclosed within a sphere of 
radius r 200m, in which the mean enclosed matter density is equal 
to 200 times the mean universal matter density p m (z) at the clus¬ 
ter’s redshift, and (2) critical overdensity masses M^ooc ( M 200 C ), 
which are analogous to (1) but enclosed within rsooc (t^ocic), and 
defined with respect to the critical density p a it(z). The correspond¬ 


ing overdensities are A 200111 , Asooc and A 2000 The critical den¬ 
sity is pcrit(a) = 3H 2 (z)/8nG, where H(z) is the Hubble pa¬ 
rameter. The mean matter density is pm(z) = Qm{z)p a «{z) with 
O m (z) = fim(l + z) 3 /E 2 (z), and where E{z) = H(z)/Hq. 


2 SIMULATIONS AND CLUSTER SELECTION 

We will refer to our hydrodynamical simulations and HMFs as 
“Hydro”, and to their dark matter only counterparts as “DMonly”. 


2.1 The Magneticum simulations 

In this work, we use a subset of the cosmological Hydro and 
DMonly boxes from the Magneticum simulation set (Dolag et al., 
in prep.) as highlighted in Table |T| The simulations are based on 
the parallel cosmological TreePM-SPH code P-GADGET3 (jSpringel 


2005). We use an entropy-conserving formulation of SPH (Springel 


& Hernquist 2002| and a higher order kernel based on the bias- 


corrected, sixth-order Wendland kernel ( Dehnen & Aly 2012} with 
295 neighbours, which together with a low-viscosity SPH scheme 
allows us to properly track turbulence within galaxy clusters l [Dolag| 
|et al.|2005||Donnert et al.|2013HBeck et al.|2016| (. 

We also allow for isotropic thermal conduction with 1/20 
of the classical Spitzer value ( |Dolag et al.|2004) . The simulation 
code includes a treatment of radiative cooling, heating from a uni¬ 
form, time-dependent ultraviolet background and star formation 
with the associated feedback processes. The latter is based on a 
sub-resolution model for the multiphase structure of the interstellar 
medium {Springel & Hemquist|2003) . 

We compute radiative cooling rates following the same pro¬ 
cedure as presented by |Wiersma, Schaye & Smith| {2009) , and 
account for the presence of an evolving ultraviolet background 
{Haard t & Madau 2001). Contributions to cooling from each 
element have been pre-computed using the publicly available 
CLOUDY photo-ionisation code (Ferland et al. 1998) for an opti¬ 
cally thin gas in (photo-)ionisation equilibrium. 

Our simulations also incorporate a detailed treatment of stel¬ 


lar evolution and chemical enrichment following Tornatore et al. 


12007 1 , a multiphase model for star-formation (Springel & Hern 


|quist 2003jl, and feedback processes associated with supemovae 
driven galactic winds and AGN {Springel & Hemquist||2003| |Di| 
|Matteo et al.|2008l|Fabian|20 1 0) . Additional details about the sim¬ 
ulation code are available elsewhere (e.g. |Hirschmann et al.|2014| . 

Initial conditions are created at redshift z m \ = 60 applying the 
Zel’dovich approximatioij^] for a spatially flat ACDM cosmology 
with parameter values from WMAP7 {Komatsu et al,|20ll) : matter 
density fl m = 0.272, baryon density fib = 0.0456, scalar spec¬ 
tral index n s = 0.963, variance in the matter fielcj^Jag = 0.809, 
and Hubble constant Ho = 70.4 km s -1 Mpc” 1 . The simulation 
boxes have periodic boundary conditions and are initially occupied 
by an equal number of gas and dark matter particles. Their rela¬ 
tive masses reflect the global baryon fraction flb/flm. To minimize 
numerical differences between the Hydro and the DMonly simula¬ 
tions, we set up the DMonly simulations with an equal number of 
two types of collisionless particles, whose masses are equal to the 
mass of the dark matter and of the gas particles, respectively, of the 
corresponding Hydro runs. 


1 We discuss the impact of the initial conditions in Appendix|B| 

2 See Equation|2]for the definition. 
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The halo mass function 3 


Table 1. Boxes of the Magneticum simulations used in this work. The number of haloes _/V( 2 _ 0 ) refer to the Hydro runs, and M200m- 


Box 

Size 7/box 

grav. 

softening 

length (kpc) 

-^particles 

^DMparticle 

-^Thalo, min 

Simulation 

N 

II 

O 



DM 

gas 

stars 


(Mq) 

(Mq) 

Hydro 

DMonly 


4/uhr 

68.1 Mpc 

1.4 

1.4 

0.7 

2 x 576 3 

5.3 x 10 7 

6.2 x 10 11 

/ 

/ 

-0 

966 

3/hr 

182 Mpc 

3.75 

3.75 

2 

2 x 576 3 

9.8 x 10 8 

1.1 x 10 13 

/ 

/ 

2/hr 

500 Mpc 

3.75 

3.75 

2 

2 x 1584 3 

9.8 x 10 8 

1.1 x 10 13 

/ 


21528 

2b/hr 

909 Mpc 

3.75 

3.75 

2 

2 x 2880 3 

9.8 x 10 8 

1.1 x 10 13 

/ 


M 

1/mr 

1274 Mpc 

10 

10 

5 

2 x 1526 3 

1.9 x 10 10 

2.2 x 10 14 

/ 

/ 

8528 

0/mr 

3818 Mpc 

10 

10 

5 

2 x 4563 3 

1.9 x 10 10 

2.2 x 10 14 

/ 


227400 


a The Hydro simulation of Box4/uhr is only run to z = 0.13. Its DMonly counterpart contains 835 haloes at redshift z = 0. 
^ Box2b/hr is only run to z = 0.2. At this redshift, it contains 109 578 haloes. 


• 

Box 1/2 z = 0.07 

4 

Gonzales 2013 

• 

Box 0/2b z = 0.90 

+ 

Kravtsov 2014 

4 

Lagana 2011 

+ 

Chiu 2014 «z) =0.9) 



M5 00 c/Mq 


Figure 1. Baryon fractions for a subset of haloes extracted from our simula¬ 
tions and for observations at low and high redshift. As discussed in the text, 
the offset between the low- and high-redshift observations is largely due to 
systematic differences in the cluster binding mass measurements. 


The Magneticum Hydro simulations have been shown to cor¬ 
rectly reproduce AGN luminosity functions jHirschmamn^et^jir 


2013 


20141 and cluster pressure profiles ( Planck Collaboration et al. 
|McDonald et al.|[2014). In Fig. |T| we compare the baryon 


fraction /bar of haloes extracted from our simulations with recent 
observational results on galaxy cluster and group mass scales. Sim¬ 
ulation values at z = 0.9 (red) are higher than the simulation values 
at 2 = 0.07 (blue), indicating that the depletion of the intracluster 
medium is increasing over time. Also shown are observational re¬ 
sults at low redshift ( |Lagana et al.||201 F[ |Gonzalez et al.||20l3] >, 
where we update the total and stellar masses for a subset of clus¬ 
ters as presented in |Kravtsov, Vikhlinin & Meshscheryakov]{20T4| , 
and at high redshift ( 2 ) = 0.9 (Chiu et al. 20161. The low red¬ 
shift observations and simulations are in good agreement, but the 
simulations exhibit systematically larger baryon fractions than are 
observed for the sample of clusters at higher redshifts. 

The comparison of the observations with the simulations is 
complicated, given that the observations rely on measured cluster 
masses. As noted in |Bocquet et al.| ( |2015) >, the X-ray hydrostatic 
masses tend to be lower in comparison to masses obtained using 
velocity dispersions or through abundance matching that includes 
constraints from external cosmological probes; other studies using 
weak lensing mass constraints come to similar conclusions (e.g. 
|von der Linden et al.|2014[|Hoekstra et al.|2015[>. In Fig. [T]the low 


redshift baryon fractions rely on X-ray hydrostatic masses, and the 
high redshift measurements adopt the abundance matching masses 
referred to above. In moving from hydrostatic mass to abundance 
matching mass, there is an increase of ~44% in the binding mass 
and a reduction of ~27% in the baryon fraction ( |Chiu et al.|2016| . 
The observed and simulated baryon fractions are in agreement if 
we adopt these shifts as the current scale of the systematic uncer¬ 
tainties. 

2.2 Halo selection 

Haloes are initially identified through a parallel FoF algorithm with 
linking length b = 0.16. The FoF links over dark matter particles 
only. We then compute SO masses (for A 2 oom, A 200 C and Asooc) of 
each halo centered at the deepest potential point with the parallel 
SUBFIND algorithm ( |Springel et al,|[2001| |Dolag et al.||2009] l. In 
Appendix [C] we discuss an alternative approach where subhaloes 
are identified as individual objects. 

To ensure that haloes extracted from the Hydro simulations are 
not affected by issues related to resolution and numerical artefacts, 
we apply very conservative convergence criteria. For each box, and 
for each overdensity A, we only consider haloes that contain more 
than 10 4 dark matter particles within ta- We then construct cata¬ 
logues applying the lower mass limits shown in Table]!] We further 
apply an upper mass limit that corresponds to the lower limit of 
the next larger box, or to 10 16 Mq for the largest boxes (see also 
Fig. i We extract cluster catalogues at seven redshift^ that are 
roughly equally spaced in cosmic time with At ~ 1.6 Gyr. This 
time step is chosen to be larger than the typical dynamic time of a 
halo, and we therefore work under the assumption that there is no 
correlation between the different snapshots. 


3 ANALYSIS METHOD 

We provide the theoretical background on the HMF and introduce 
the fitting form we adopt. We also present the method used to per¬ 
form the multi-dimensional fits for the HMF parameters. 


3.1 The halo mass function 


The comoving number density of haloes of mass M is 

dn _ , . p m din a -1 

dM ~ I(a) M dM ’ 


0 ) 


3 We use the Hydro run of Box4/uhr for redshifts 2: > 0.13. 
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with the mean matter density p m (at redshift a = 0), and 

a 2 (M,z) = J P{k,z)W 2 (kR)k 2 dk, (2) 

which is the variance of the matter density field P(k, z) smoothed 
with the Fourier transform W of the real-space top-hat window 
function of radius R = (3M/47rp m ) 1 ' 3 ■ The function f(a) is com¬ 
monly parametrized as 


f(v) = A 


+ 1 


exp 




(3) 


with four parameters A, a, b, c that need to be calibrated (e.g. |Jenk-| 
|ins et al.|2 001 ). Flere, A sets the overall normalization, a and b are 
the slope and normalization of the low-mass power law, and c sets 
the scale of a high-mass exponential cutoff. 

The fit parameters in f(a) will in general depend on redshift 
and cosmology, which must be accounted for. In the following we 
note that this dependence is weak in certain special cases; this is 
referred to as universality of the HMF. 


3.2 Halo mass function for spherical overdensity masses 

Many studies of the HMF are performed using the FoF technique. 
For a linking length b ~ 0.2, the resulting HMF is very close to 
being universal ( [Jenkins et al.|200l) . However, for observational 
reasons, real cluster masses are measured in terms of SO masses. 


3 . 2.1 The halo mass function for A200 m 

When using a suitable SO A mea n ~ 180, the above HMF fitting 
function is still close to being universal ( [Jenkins et al.|200l) . Sim¬ 
ilarly, Tinker08 use A200111 as their universal mass definition, and 
Watsonl3 argue for Ai7g m . These overdensity definitions are all 
very similar; we adopt A200111 in this work. 

We follow previous work and allow for departures from uni¬ 
versality by parametrizing a possible redshift dependence of the 
HMF in A200111 as a power law in 1 + 2: 

A(z) = -A 0 (l + z) A * 
a(z) = a 0 ( 1 + z) a * 
b(z ) = 6o(l + z) bz 

c(z) = co(l + z) Cz (4) 

where the subscript 0 denotes the values at redshift 2 = 0, and 
where A z ,a z ,b z , c z are additional fit parameters. Note that many 
authors assume the cutoff scale c to be constant (c z = 0, e.g. Tin- 
ker08, Watsonl3). 


3.2.2 Toward a universal mass function for Asooc and A200C 

We also wish to calibrate the HMF for Asooc, which is a convenient 
mass definition within X-ray studies of clusters where the emis¬ 
sion cannot easily be traced beyond rsooc, and for A 200 C, which is 
used for measurements of cluster galaxy velocity dispersions and 
for weak gravitational lensing shear profiles. It is not a priori clear 
that one can simply use the same form of the fitting function that is 
valid for A200111, and calibrate it against simulations at other over¬ 
densities, as one might miss some redshift and cosmology depen¬ 
dent behaviour. Remember, for example, the very different redshift 
evolution of p m (z) and p cr it(z). 

Tinker08 provide the HMF for a range of different A mean . One 
then uses these HMFs together with A mea n(z) = A cr u/fl m (z) to 


convert from critical to mean density as a function of redshift. Their 
approach relies on the implicit assumption that the fitting function 
correctly captures the behaviour for every A mean . Watson 13 provide 
a correction to their Ai 78 m HMF that depends on A mean (z). 

We choose a novel approach, based on the assumption that 
the HMF is approximately universal in A 2 oom- Using the HMF for 
A 200111 we then determine the HMF for Asooc as: 
dn _ dn dM2oom 
dM^ooc rfM200m dM$ooc 

p m rflng -1 \ M 5 oq c 
AT50OC dMsoOc J MiOOm 



The first factor now has the same functional form as the HMF for 
A200111 (Equations [TJQ. The second factor describes the conver¬ 
sion between the two different SO mass definitions. This conver¬ 
sion is discussed below. The crucial point is that the new functional 
form of the HMF established in Equation [5] has the same close-to- 
universal properties as the HMF for A 2 oom- Note that, due to un¬ 
certainties and systematic scatter in the ATsooc-A^oom conversion, 
/(a) in Equation [ 5 ] still needs to be calibrated against the Asooc 
simulation data. 

SO masses can be converted from one to the other assuming 
a cluster density profile and a mass-concentration relation. We use 
the parametrizations by |Navarro, Frenk & White~| < [ 1997[ ) and |Duffy| 
et al. |2008| and establish an analytic fit for ATsooc/A^oom- This 


fit for A/ 500 C/AT 200111 depends on mass, redshift, and fl m (which is 
involved in the overdensity conversion). The following prescription 
reproduces the mass conversion at the few percent level in the range 
0 < 2 < 2, 10 13 < Msqoc/Mq < 10 16 , and 0.1 < O m < 0.5: 


Msooe 

M 2 00m 


a + /Jin 


Msooe 


The parameters a and /3 are functions of fl m and redshift: 

/ 3 (H m ) = -1.70 x 1CT 2 + H m 3.74 x 1CT 3 


z) 


at z + a 2 

a 0 -;- 

z + a 2 


a 0 (H m ) = 0.880 + 0.329 
ar(n m ) = 1.00 + 4.31 x 10" 2 /12 m 
a 2 (fi m ) = -0.365 + 0.254/n m . 


(6) 


(7) 


In an analogous way, we establish the form of the HMF for 
A 200 C in Appendix[A] In Section[6]we present step-by-step instruc¬ 
tions of how to compute our HMF for either of the three SO defini¬ 
tions. This presentation explicitly shows how the SO mass conver¬ 
sions are used. 


3.3 Finite volume correction 

Throughout this work, we use cluster samples produced by simula¬ 
tions to understand the HMF observed in the real Universe. How¬ 
ever, there is one subtle difference that needs to be accounted for: 
in contrast to the Universe, every simulation box is finite in size. 
Therefore, we can only capture modes in the density field that are 
smaller than Lbox- This means that there is an upper mass limit 
corresponding to the longest modes, beyond which the simulations 
will systematically underestimate the number of objects. 

We correct for this effect following the approach of previous 
analyses ( |Lukic et al,|2007[|Bhattacharya et al.|2011| Watsonl3). 
Briefly, the variance of fluctuations a(M) is corrected by subtract¬ 
ing the variance at scales corresponding to the box size a(Rbox)'- 

(^corrected (AT) = aleory(M) - a 2 (R box ), (8) 
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Figure 2. Halo mass function dN/dlnM from our Hydro and DMonly simulations. Redshift is increasing from top to bottom and takes values z = 
0, 0.13,0.3, 0.5, 0.8.1.2, 2. The data points are slightly offset in mass for improved readability. We also show the 2a allowed regions of our Hydro (blue 
region) and DMonly (green region) fits. For Hydro, the error bars and allowed region are narrower because we also use data from the much larger BoxO/mr. 


where, for simplicity, we equate the spherical volume 4/3-7rft^ ox to 
the cubical simulation volume L^ ox . However, as we apply upper 
mass limits to the cluster samples (see Section [2~2] i, the correction 
has negligible impact on our analysis. In fact, for each box size, 
the correction would become important at masses that are about 2 
orders of magnitude larger than the corresponding upper mass limit 
we apply. We test these finite volume effects by reproducing the 
results presented in Section [4] without the correction; in this case, 
the results do not significantly change. Nevertheless, we apply the 
correction to each of our boxes. 


3.4 Parameter estimation 

We use a Bayesian likelihood approach, which allows us to cor¬ 
rectly capture the Poisson errors on the measured number of clus¬ 
ters as a function of their mass and redshift. This choice differs from 
using (Gaussian) \ 2 statistics and jackknife errors (e.g. Tinker08), 
or corrections to x 2 statistics to account for the Poisson errors (e.g. 
Watson 13). 

The likelihood at each point p in parameter space is calculated 
in the following way: We calculate the matter power spectrum using 
the transfer function of Eisenstein & Hu (1998 1999), taking bary- 
onic effects into account. This is the same prescription used to set 
up the initial conditions of the Mcigneticwn simulations. We eval¬ 
uate the likelihood C by applying Poisson statistics in log-spaced 


mass bins of size A log 10 M = 0.1 ( jCas h| 1979j ): 

, r _ ST' i dn(Mi\p,Zi) f dn(M\p, z) 

n ^ dM J dM 


dM, (9) 


up to an arbitrary constant, and where i runs over all clusters in 
the sample. The second term equals the total number of expected 
clusters. We have checked that decreasing the bin size does not 
change our results. 

In practice, given a set of parameters p, we perform the above 
calculation for each redshift and for each of the simulation boxes, 
and sum the log-likelihoods. When fitting for the HMF in this way, 
we are facing a problem with moderately large dimensionality (8 
free parameters); we utilize the emceaHcode for efficient explo¬ 
ration of parameter space ^Foreman -Mackey et al.|20l3) . We test 
our fitting procedure against several mock catalogues that contain 
a factor 100 times more clusters than our simulation data. In these 
tests we recover the input values within the statistical uncertainties 
and conclude that our fitting method is unbiased to a level that is 
much smaller than the uncertainties we report. 


4 MAGNETICUM MASS FUNCTIONS 

In Fig. [2] we show the number density of haloes from our Hydro 
and DMonly simulations at seven different redshifts. The error bars 

4 http://dan.iel.fm/emcee/current 
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Figure 3. Number density of haloes relative to our Hydro fit for M 200 m (left panels), M 200 c (center panels), and M 500 C (right panels). The coloured bands 
correspond to the 2cr allowed regions of our fits; the data points are slightly offset in mass for better readability. For masses above M > 2.2 X 10 14 Mq, the 
Hydro (DMonly) data points shown here are extracted from BoxO (Boxl). The red dashed line shows the fit by |Watson et al.|(20 1 3| , and the orange line shows 
the fit by |Tinker et al.1(2008) . 


show the Poisson uncertainty on the measured numbers. We also 
show the 2cr allowed regions of the fitting functions. The charac¬ 
teristic low-mass power law and the steepening at high mass are 
clearly visible. 

The impact of baryons can be better distinguished in Fig. [3] 
which shows the number density of clusters from both our simula¬ 
tions relative to the best-fitting function for Hydro. We also show 


the allowed 2a regions of our fits and two external dark matter only 
fitting functions (Tinker08, Watson 13). 


4.1 Impact of baryons 

Throughout the different SO definitions, the comparison of our Hy¬ 
dro and DMonly simulations tells an interesting story: Atz = 2, the 
highest redshift we consider, there is essentially no difference be- 
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The halo mass function 7 


Table 2. Best-fitting HMF parameters for DMonly (Box 1+3+4) and Hydro (BoxO+1+3+4) simulations. The covariance matrix for Hydro M 5 00c is shown in 
Table[3] The HMFs for M 500 C and M 200 C are calculated from Equations |5land|AT| respectively. 


Parameter 

A 

a 

b 

c 

A z 

a z 

b z 

C z 

AT2oom 

DMonly 

0.175 

1.53 

2.55 

1.19 

-0.012 

-0.040 

-0.194 

- 0.021 

Hydro 

0.228 

2.15 

1.69 

1.30 

0.285 

-0.058 

-0.366 

-0.045 

AT 200 C 

DMonly 

0.222 

1.71 

2.24 

1.46 

0.269 

0.321 

-0.621 

-0.153 

Hydro 

0.202 

2.21 

2.00 

1.57 

1.147 

0.375 

-1.074 

-0.196 

ATsooc 

DMonly 

0.241 

2.18 

2.35 

2.02 

0.370 

0.251 

-0.698 

-0.310 

Hydro 

0.180 

2.29 

2.44 

1.97 

1.088 

0.150 

-1.008 

-0.322 


Table 3. Covariance matrix for our Hydro HMF in M 5000 The corresponding best-fitting values are 

A a b c A z 

shown in Table [ 2 ] 

a z 

b z 


Cz 

A 

3.33 x 10~ 4 

4.58 x 10" 4 

-2.61 x 

10~ 3 

1.89 x 

lO" 5 

-3.56 x 

10 -3 

-2.25 x 

10~ 4 

2.22 x 

10 -3 

7.38 x 10" 5 

a 


6.44 x 10" 3 

-7.72 x 

10~ 3 

1.60 x 

10” 3 

3.31 X 

10“ 3 

2.67 x 

10“ 5 

-3.55 x 

lO" 4 

-4.13 x 10“ 4 

b 



2.39 x 

10 -2 

-1.10 X 

10“ 3 

2.34 x 

10“ 2 

1.88 X 

10“ 3 

-1.63 x 

10“ 2 

-3.66 x 10" 4 

c 





5.12 x 

lO" 4 

2.71 x 

10“ 3 

3.01 x 

lO" 4 

-1.46 x 

10“ 3 

-2.05 x 10“ 4 

A z 







1.16 x 

lO" 1 

1.65 x 

10" 2 

-7.71 x 

10 -2 

-3.08 x 10" 3 

a z 









4.18 x 

lO" 3 

-1.24 x 

10 -2 

-2.62 x 10“ 4 

b z 











5.33 x 

10“ 2 

2.10 x 10“ 3 

C z 













2.40 x 10“ 4 


tween our Hydro and DMonly HMFs. As structure formation con¬ 
tinues, baryonic effects become important, and at redshift z ~ 1, 
we observe that baryon depletion is important for low-mass clusters 
up to ~ 10 14 Mq. At fixed mass, this reduces the number density 
of clusters by < 10%. Further following the redshift evolution, the 
effects of baryon depletion propagate up to more massive clusters. 
For low redshifts z < 0.3, and for masses up to ~ 10 14 AfQ, our 
Hydro HMF is low by about 10 — 30% compared to the DMonly 
case, and the difference increases when going from A200111 to A200C 
to A 5 ooc- This dependence on the SO definition is consistent with 
the picture of baryons being ejected from the cluster central re¬ 
gions, as r 50 0c < f200c < t*200m- 

In the mass range between 10 14 — 10 15 Mq, there is essen¬ 
tially no difference between our Hydro and DMonly simulations; 
as expected, baryonic effects are negligible on these scales. For 
even higher masses > 10 15 Mq, and for A200111 and A200C, the 
Hydro simulations contain fewer clusters than DMonly. First, we 
check that our simulations have converged and correctly estimate 
the abundance of the most massive haloes. As discussed in Ap¬ 
pendix [B] we run two additional DMonly realizations of Boxl/mr, 
initialized at higher redshifts zm = 80 and z\m = 100, and find that 
they agree very well with our default simulations initialized at red¬ 
shift Zi n i = 60. Then, the high-mass end is much better constrained 
for Hydro where we also use data from BoxO/mr which is not avail¬ 
able for DMonly. We therefore expect that these differences would 
reduce when analyzing a larger simulation box for DMonly. 


4.2 Comparison with other dark matter only fits 

The comparison of our simulations and fits with Tinker08 and Wat¬ 
son 13 can be summarized as follows: 

(i) At A 200111 , the fit by Tinker08 and our DMonly fit are in very 
good agreement. This is an important cross-check, since both simu¬ 


lations assume the same physical model and are similar in volume, 
and the assumed functional forms of the HMF are identical. 

(ii) The fit by Watson 13 tends to predict less low-mass haloes 
than the other HMFs considered here. At the high-mass end, for 
A 200111 , Watson 13 predict less haloes than Tinker08 or our fits, but, 
for larger SO, Watsonl3 roughly agrees with Tinker08. 

(iii) For the larger SO A 200 C and Asooc, Tinker08 predicts up 
to ~ 25% more haloes for AT < 1O 12 M 0 ; this difference steadily 
reduces to zero for masses M > 10 14 ATq. 


Item (iii) above is due to a different approach in halo find¬ 
ing: We identify haloes and extract their masses using SUBFIND. 
The SO masses are constructed around the point corresponding to 
the minimum in the gravitational potential, neglecting substructure 
that eventually exceeds the overdensity A when going to high val¬ 
ues like Asooc. In contrast, Tinker08 use an SO finder that reveals 
these subhaloes. In Appendix [C] we confirm that the presence of 
substructure increases the abundance of objects by about 25% on 
low mass scales M ~ 10 12 Mq\ the effec t drops to zero for larger 
masses Msooc 10 13 ATq (see also Kravtsov et al.|2004 1 . There¬ 
fore, we expect this effect to be negligible for current cluster stud- 


4.3 Halo mass function fits 

In Table [2] we present the best fit parameters for the fits to the 
data from our DMonly and Hydro simulations. We also evaluate 
the x 2 P er degree of freedom v at the best-fitting location and the 
corresponding probability to exceed. The latter takes values that are 
very close to unity for all six HMF fits, indicating that the adopted 
form of the HMF is suitable for fitting the data. 

The HMF parameters for A 2 oom can be directly compared 
with the literature because in this case the functional form of the 
HMF is identical. For DMonly, we find a value of the exponential 
cut-off scale c= 1.19±0.03 that is fully consistent with Tinker08 
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Figure 4. Simulation of idealized representations of the SPT. Planck and 
eROSITA cluster catalogues. Top: Selection functions used to create the cat¬ 
alogues. Bottom: Distribution of simulated clusters as a function of redshift. 


(c = 1.19), but significantly smaller than the c = 1.318 in Wat- 
sonl3. Note that a large value of c corresponds to a low cut-off 
scale in mass (see Equation [3]>, resulting in the underprediction of 
massive haloes by Watsonl3 as shown in Fig. [3] Finally, there is no 
evidence for redshift-evolution of the cut-off scale in the DMonly 
case with c z = —0.021 ± 0.054, but clear evidence in the case of 
the Hydro HMF, where c z = —0.045 ± 0.012. 

The HMFs for A 200 C and Asooc are often used in observa¬ 
tional studies. For the latter SO, we also show the covariance matrix 
for the Hydro HMF parameters in Table[3] This information should 
be used with the best-fitting parameters from Table [2]to capture the 
(statistical) uncertainties related to our HMF. For reference, the la- 
uncertainty on the HMF at Msooc = 10 15 Mq is ~ 1%. 


5 COSMOLOGICAL IMPACT 

There are differences between the HMFs extracted from our Hydro 
and DMonly simulations, and those from the literature. When used 
to interpret real cluster samples, the different HMFs will ultimately 
lead to different cosmological results. In the following, we quantify 
and discuss this effect. To this end, we create simulated cluster cata¬ 
logues using our best-fitting Hydro HMF, and use either the Hydro, 
the DMonly, or literature fits to perform cosmological analyses. Be¬ 
cause the baryonic impact on the HMF depends on mass and red- 
shift, we expect qualitatively different shifts when using different 
HMFs, depending on the properties of a specific cluster survey. 

We create and analyse three sets of simulated catalogues, 
whose properties approximately match real samples from the South 
Pole Telescope (SPT; [Carlstrom e t al.||2011|>, the Planck satellite 
( |Tauber et al.|20 1 0~} and eROSITA l |Predehl et al.|2014) . The selec¬ 
tion functions we assume are shown in Fig.QJand will be discussed 
in more detail. All samples are defined for SO A 5000 

We use the fit method described in Section[3] but we now fit for 
the cosmological parameters, and keep the HMF parameters fixed. 


Since no covariance matrix is available for the literature HMFs we 
compare to, we use our best-fit parameters without uncertainties, 
too, in order to make a comparison on equal footing. We further 
show that the (statistical) uncertainty on our Hydro HMF has negli¬ 
gible impact on the cosmological constraints. We restrict this analy¬ 
sis to the parameters D ra and og, which strongly affect the measured 
cluster abundance. Remember that these parameters enter the HMF 
calculation in Equation [I] through their impact on the matter power 
spectrum P(k, z) and the matter density p m . The D m -og likelihood 
contours from the cluster number counts experiment exhibit a char¬ 
acteristic, elongated degeneracy in the D m -o-g plane (see Figs[5]and 
[6|. The parameter combination o-g(D m /0.27)°' 3 is interesting be¬ 
cause it reflects the width of this degeneracy, i.e. the direction in 
D m -o- 8 space which is best constrained using clusters. We show the 
constraints we recover on this parameter combination, too. 

In this test, we directly use the simulated cluster masses. That 
is, we do not include any systematic uncertainties and measurement 
errors related to mass estimation as one would have to do for a real 
cluster sample. This also means that the uncertainties we recover 
only represent the statistical uncertainties related to the sample size, 
and cannot be compared with results from observed clusters. The 
aim of this analysis is to investigate and quantify offsets related to 
the HMF, which justifies this simplified approach. For this same 
reason, we do not quote the errors on the recovered parameters. 

The typical uncertainties on the cosmological parame¬ 
ters from current cluster samples_are q-( D m ) ~ a(aA 
a(a 8 (D m /0.27)°' 3 ) ~ 0.03 (e.g. 


Planck Collaboration et al. 


2014 


|Mantz et al.|20 1 5) . We will refer to these characteristic numbers ii 
the following. 


5.1 Cosmological analysis of an SPT-like cluster sample 


The SPT sample is selected through the cluster Sunyaev- 
ZeTdovich effect (SZE; |Sunyaev & ZeTdovich||l972] > signature, 
and we approximate the catalogue as mass-selected with Msooc > 
3 x 10 14 M(7), an d restrict to redshifts z > 0.25 (see Fig.[4]and cf. 


Bleem et al. 


2015). For the SPT survey of size 2500 deg - ', our sim¬ 


ulated catalogue contains 687 systems. We consider three different 
input cosmologies with different values of D m = 0.22, 0.272, 0.4 
with the same as = 0.809 in each case. A subset of the results 
appears in Tableland Fig.[5] 

The results from both our Hydro and DMonly HMFs show 
nearly perfect agreement, indicating that the effect of baryons on 
the HMF is negligible in this case. This is expected, because, as 
previously noted, the impact of baryons is most important for low- 
mass clusters at low redshifts (see Section[4Tj. The SPT-like sur¬ 
vey does not probe this mass and redshift regime. 

The marginalized constraints on a 8 obtained using Tinker08 
are in very good agreement with the ones from our HMFs, but 
Dm slightly shifts by —0.01. Due to a mild tilt in degeneracy axes, 
o-g(D m /0.27) 0,3 is slightly tighter when using the Tinker08 fit. 

The constraints obtained on D m using the Watson 13 fit are 
tighter than the ones just discussed (see Fig. [5}. However, these 
results seem to be biased toward D m ~ 0.27. For example, the 
preferred value recovered for the sample with input D m = 0.22 is 
D m (Watsonl3) = 0.25, and we further obtain D m (Watsonl3) = 
0.31 for an input value D m = 0.4. Their assumed form of the red¬ 
shift dependence of the fit parameters (Equation [4] in this work. 
Equations 13-15 in Watson 13) involves D m (s). We suspect that this 
parametrization introduces an implicit and spurious preference for 
D m ~ 0.27, which is their simulation input value. We will not con¬ 
sider the Watsonl3 fit for the rest of this work. 
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(a) Results for input Q m = 0.272. 


(b) Results for input Q m = 0.22. 


Figure 5. Cosmological results from simulated SPT-like cluster samples, created using our Hydro halo mass function. The two different input cosmologies are 
marked by the dashed lines, in both cases erg = 0.809. We show likelihood contours (68% and 95% confidence) in f2 m -(rg-(Tg(Q m /0.27) 0 ' 3 space. Results 
obtained from different halo mass functions are colour-coded. The baryonic impact of the halo mass function is negligible for this high-mass and h igh-re dshift 
sample. For the input value fl m = 0.272, there is a small shift toward smaller values of D m and <rg(S2 m /0.27)°' 3 when using Tinker et ak (2008 1 . The 
constraints obtained from|Watson et al.~|(20l 3| are tighter, but seem to be biased towards f! m ~ 0.27. See discussion in the text for more details. 


Table 4. Mean recovered cosmological parameter values from different halo 
mass functions for simulated cluster samples. We do not show errors be¬ 
cause we are interested in absolute shifts only. 


Parameter 


r^m 

O- 8 

(T8(n m /0.27) 0 - 3 

Input 


0.272 

0.809 

(0.811) 

SPT-like sample 




Hydro 


0.272 

0.808 

0.810 

DMonlv 


0.269 

0.809 

0.808 


Tinker et al. 

|2008 

0.260 

0.813 

0.804 

eROSITA-like sample 



Hydro 


0.272 

0.809 

0.811 

DMonlv 


0.260 

0.812 

0.803 


Tinker et al. 

|2008 

0.258 

0.815 

0.804 

Planck -like sample 



Input 


0.316 

0.830 

(0.870) 

Hydro 


0.320 

0.828 

0.871 

DMonlv 


0.312 

0.836 

0.873 


Tinker et al. 

|2008 

0.314 

0.816 

0.853 


The statistical uncertainty of our Hydro HMF is captured by 
the covariance matrix in Table [3] We repeat the cosmological anal¬ 
ysis using the Hydro HMF and its parameter covariances, and in¬ 
fer the additional uncertainties due to the uncertainty on the HMF 
(thmf using quadrature addition. We find <thmf(Hhi) ~ 0.004 and 
ohmfOs) ~ 0.002, and conclude that the statistical uncertainties 
on our HMF are completely negligible for current cluster samples. 


5.2 Cosmological analysis of a Planck -like cluster sample 

The Planck cluster sample is selected using the SZE, too, and ex¬ 
tends down to redshift 2 = 0. However, the satellite’s beam is larger 
than the SPT beam, and the survey mass limit varies significantly 
with redshift. We mimic the Planck selection function following the 
sample mass-redshift distribution shown in Fig. 1 in |Planck Col-| 
|laboration et al.| j2015b| >. We further assume a hydrostatic bias of 
1 — b ~ 0.8, which then leads to the selection function we show in 
Fig. 0 For this exercise, we choose our input cosmology to match 
the values preferred by the Planck CMB anisotropy measurement 
(H m = 0.316, erg = 0.83, H 0 = 67.3 km s’ 1 Mpc' 1 ; 
[Collaboration et al.|2QI5afr . Assuming a sky coverage of 65%, the 
simulated catalogue contains 639 clusters. 

The results are shown in Table [4] and Fig. |6(a)| We recover 
very similar constraints on S2 m for our Hydro and DMonly HMFs 
and Tinker08. There is some spread in the constraints on as, where 
the difference between our DMonly and Tinker08 is A(crg) ~ 
0.02. We recover identical constraints on crg(H m /0.27)°' 3 from 
our Hydro and DMonly fits, but using Tinker08 leads to an shift 
A(crg(f2 m /0.27)°' 3 ) ~ —0.02. The shift in ag and the induced off¬ 
set in crg(n m /0.27)°' 3 are larger than the differences in our analy¬ 
sis of the SPT-like sample. As discussed in Section |4~2| our HMFs 
and Tinker08 differ somewhat at the highest clusters masses for 
redshifts 2 > 0.3. Then, as shown by the selection function in 
Fig. 0 the Planck-sample contains more massive clusters than the 
SPT-sample which makes the Planck -sample more sensitive to dif¬ 
ferences in the high-mass end of the HMF. 

The analysis of the real Planck cluster sample is limited 
by the systematic uncertainty on the overall mass scale, often 
parametrized by the mass bias factor. In their latest cluster cos¬ 
mology analysis, the Planck collaboration adopts measurements of 


Planck 


© 0000 RAS, MNRAS 000, 000-000 










































10 S. Bocquet, A. Saw, K. Dolag and J. J. Mohr 




a 8 (Q m /0.27) 0 - 3 


(a) Results for a simulated Planck -like cluster sample. The input cosmol¬ 
ogy is chosen to match the Planck CMB constraints {Planck Collaboration! 
|et al.|2015a| . Given the large mass limit of the sample there is little dif¬ 
ference between the results from our Hydro or DMonly halo mass func¬ 
tions. Using the[Tinker et al.||2008j> fit leads to shifts A(crs) ~ —0.01, and 
A (erg (n m /0.27) 0il ) ~ —0.02. 


(b) Result for a simulated eROSITA-like cluster sample. For this sample, 
there is a difference between the results from the Hydro and the DMonly 
halo mass functions, which is due to the impact of baryons. Neglecting the 
baryonic impact on the halo mass function leads to an underestimate of 
A(Q m ) ~ A(tr g (n m /0.27)°' 3 ) ~ -0.01. The DMonly and the|Tinker| 
|et al.|(2008j fits produce very similar constraints for this sample. 


Figure 6. Cosmological results from simulated realizations of approximations to the Planck and eROSITA cluster samples, created using our Hydro halo mass 
function. The input cosmologies are marked by the dashed lines. We show likelihood contours (68% and 95% confidence) in Q m -crg-crg(Q m /0.27) 0 ' 3 space. 
Results obtained from different halo mass functions are colour-coded. Note the very different scales in both panels. 


the mass bias from different authors i j Planck Collaboration et ah] 
|2015b[ l. Assuming a measurement of the mass bias from ground- 
based weak tensing, the recovered values of erg are low by roughly 
0.05 — 0.08 compared to the value recovered from the Planck 
CMB anisotropy measurements. However, due to the combination 
of systematic and statistical uncertainties, the Planck CMB pre¬ 
ferred value in flm-os space is still compatible with the cluster 
measurement at the ~ 1 — 2<r level. Our simplified analysis sug¬ 
gests that using our HMF instead of Tinker08 could lead to a shift 
in <T8(f2 m /0.27) 0 ' 3 by about +0.02, which would reduce the dif¬ 
ferences between the values preferred by Planck clusters and the 
CMB anisotropy. 


5.3 Cosmological analysis of an eROSITA-like cluster sample 

The eROSITA cluster sample will be X-ray flux selected and ex¬ 
tend from redshift z = 0. For the present test, we assume a detec¬ 
tion limit of 50 photons in the 0.5 — 2.0 keV band with a typical 
exposure time of 1.6 ks. We model this selection as a combination 
of a redshift-dependent mass threshold Msooc >2.3 zx 10 14 M©, 
with an additional mass cut at 7 x 10 13 Mq (see Fig. Mj and com¬ 
pare with fig. 2 in|Pillepich, Porciani & Reiprich|(|2012{, and also 
|Merloni et al~| < ]2012[ >; |Borm et al.| < |2014) ). The eROSITA full-sky 
catalogue simulated in this way contains ~ 1.45 x 10 s clusters. 

The results of the analysis of this sample appear in Table [4] 
and Fig. |6(b)| The recovered constraints are very tight due to the 
large cluster sample and the fact that we do not include mass mea¬ 
surement uncertainties. For this sample, the constraints on erg agree 
to within A(crg) ~ 0.006. More importantly, the recovered con¬ 


straints on O m are in good agreement between our DMonly and 
Tinker08, but both are low compared to the results obtained us¬ 
ing the Hydro HMF. This is an indication that baryonic effects are 
indeed important for this sample. As previously discussed, baryons 
have their strongest impact on the HMF at low redshifts and for low 
masses, which is a regime that is well probed by eROSITA. There¬ 
fore, neglecting the baryonic impact for this sample and using a 
dark matter only HMF would lead to a bias of Afi m ~ —0.01. 
This bias is of the same order as the expected uncertainty from 
eROSITA <7(f2 m ) ~ 0.012 i Pillepich, Porciani & Reiprichj2012 


meaning that the impact of baryons on the HMF will have to be 
accounted for in the cosmology analysis. For the same reasons just 
discussed, (jg(fl m /0.27) 0 ' 3 is slightly underestimated when using 
either our DMonly or the Tinker08 fit. 


6 SUMMARY 


We calibrate the HMF and investigate the impact of baryons us¬ 
ing the hydrodynamic Magneticum simulations together with dark 
matter only counterparts. Our simulations and the halo selection are 
characterized by (1) a treatment of the baryonic component and of 
AGN feedback that is in good agreement with several observations 
such as baryon fractions, AGN luminosity functions i jlfirschmann 
et al.||2014| and cluster pressure profiles (Planck Collaboration 
et al.|2013[ McDonald et al.|2014| l, (2) large cosmological volumes 
probed by boxes of up to (3818 Mpc) 3 , which allow us to track 
cluster masses up to a A 1 > 10 15 Mq, and (3) a conservative halo 
selection with > 10 4 dark matter particles within ta, minimizing 
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potential biases related to numerical resolution. To avoid a different 
sampling of the initial density fluctuations, the DMonly simulations 
were run using two species of dark matter with masses correspond¬ 
ing to those of the dark matter and baryonic particles in the Hydro 
simulations. We extract SO masses A 200111 , A 200 C, and Asoo c . 

The presence of baryons tends to decrease the cluster masses, 
which - given the shape of the HMF - leads to a decrease of the 
expected number of objects for a given mass (see Figs [2] and [3}. 
The number density of haloes decreases by up to ~ 15% for low 
masses M < 1O 14 M0 and at low redshifts z < 0.5. At higher 
masses and redshifts, our Hydro and DMonly simulations agree 
very well. Qualitatively similar results have been recently presented 


in other publications (Cui, Borgani & Murante 2014, 

Cusworth 

et ah 

2014 

[Velliscig et al.pOMJVogelsberger et al.|201^ 

Schaller 

et ah 

2015 

'. Martizzi et ah (2014) find a mild increase of the HMF 


due to baryons. 


The HMF shape varies only weakly with redshift and cosmol¬ 
ogy when masses are defined either by FoF with b ~ 0.2, or for 
SO A mcan ~ 200. Therefore, for the HMF for Asooc (A 2 ooc), we 
introduce a mapping between Msooc (Azotic) and M 200111 as a func¬ 
tion of mass, redshift, and O m , and argue that this allows us to use 
the universal properties of A 200111 also for masses defined by Asooc 
(A 2 ooc). In practice, our HMF fits are used as follows: 

• Calculate a(M, z) using Equation]!] 

• Calculate f(a, z) from Equations[3]and[4] using the appropri¬ 
ate parameters from Table [2] 

• For A 200111 : The HMF is given by Equation]!] 

• For Asooc: Calculate Msooc/A-faoom from Equations [6] and [ 7 ] 
and obtain the HMF from Equation]!] 

• For A 200 C: Calculate A^ooc/ATooom from Equations ! A2| - |A4| 
and obtain the HMF using Equation |Al| 

Note that the same approach could be used to propagate the univer¬ 
sal behaviour of the Aaoom HMF to any overdensity A. 

We investigate how the differences among our Hydro, 
DMonly and previously published dark matter only HMFs affect 
cosmological results from cluster abundance measurements. To 
this end, we simulate idealized representations of the SPT, Planck, 
and eROSITA surveys, assuming simplified selection schemes as 
shown in Fig. [4] We assume perfect knowledge of cluster masses 
A/sooc, and do not account for any uncertainties or systematics re¬ 
lated to mass-observable relations. Therefore, the cosmological pa¬ 
rameter uncertainties we recover here are tighter than the actual 
constraints that would be obtained in a comprehensive analysis of 
real data. Moreover, neglecting the conversion from observable to 
mass likely removes some cosmological dependencies. However, 
this test can be used as guidance in understanding the impact of 
differences in the HMF. 

The results of these analyses can be summarized as follows 
(see also Figs[5]and[6]and Table]!}: 

• For the SPT-like and Planck- like samples, the impact of 
baryons is negligible, and we obtain identical cosmological results 
using either our Hydro or DMonly HMFs. 

• For the SPT-like sample, results obtained using the Tinker08 
fit essentially agree with results obtained using our HMFs. 

• The HMF by Watson 13 seems to bias results toward S2 m ~ 
0.27. This may be due to their parametrization of the redshift evo¬ 
lution of the HMF shape parameters using fi m (,s), which results in 
a heightened and likely artificial cosmological sensitivity. 

• For the Planck- like sample, using our Hydro HMF instead of 
Tinker08 shifts the results by A(o8(fi m /0.27)°' 3 ) ~ 0.02. This 


shift corresponds to about half the observed difference between the 
latest Planck clusters and CMB constraints; using our HMF should 
therefore lead to better agreement between the two probes. 

• The eROSITA sample extends to lower masses than the 
SPT and Planck catalogues. We observe an offset in the results 
from Hydro and DMonly, which we identify as the impact of 
baryons. Neglecting this effect leads to an underestimate Afl m ~ 
A(crg(f2 m /0.27) 0 ' 3 ) ~ —0.01, which is comparable to the ex¬ 
pected overall uncertainty. 

Part of the differences between the cosmological results recov¬ 
ered using our HMFs and using the Tinker08 could be due to differ¬ 
ent interpretations of the connection between the observed galaxy 
clusters and their representation as haloes in the simulations. Dif¬ 
ferent assumptions lead to different halo-finding approaches: We 
extract SO masses around the minimum in the gravitational poten¬ 
tial in each halo using SUBFIND. In contrast, Tinker08 employ an 
SO finder, and also include subhaloes in their fit. These different 
approaches only affect the HMF for low masses Msooc < 10 13 Mq 
and are therefore unimportant for the cosmological study of current 
cluster samples. 

More work, both on the theoretical and on the numerical as¬ 
pects of calibrating the HMF is needed to be able to fully extract 
the cosmological information from near-future cluster samples. It 
is important to better understand the cosmological dependencies of 
the fitting functions, and to construct an analytic formula whose 
universality - or indeed departure from universality - is well un¬ 
derstood. Finally, a careful comparison of cluster catalogues gener¬ 
ated from different sets of numerical simulations would be helpful 
to better understand the systematic uncertainties on the HMF. 
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Figure Al. Number of haloes extracted at redshift 2 = 0.5 and z = 0 
from simulations started at initial redshift z m \ = 60 (used in this work) 
and two control simulations started at higher initial redshifts z m [ = 80 and 
^ini — 100- The error bars show Poisson uncertainties. All three simulations 
agree within the statistical uncertainties, indicating good convergence of the 
simulations used in this work. 


Figure Bl. Impact of substructure on the DMonly HMF for A 5000 We 
show the abundance of haloes, including subhaloes, relative to the num¬ 
ber of main haloes as used throughout this work. For small halo masses, 
substructure has a significant contribution to the abundance, and the abun¬ 
dance of haloes including substructure seems to be in good agreement with 
|Tinker et al.| {2008}. The impact of subhaloes vanishes for haloes with 
A/500C ^ 10 1 Mq. 


APPENDIX A: SPHERICAL OVERDENSITY A 20 oc 
Applying the method described in Section[L2| the HMF for M 2 00 c 


dn _ p m d\na 1 M 2 qoc 

dM 2 ooc M 2 ooc dM 2 ooc M 2 oom 


(Al) 


Assuming a cluster density profile ( |Nayarro^Frenk & White] 1997} 
and a mass-concentration relation ( [Duffy et al.|2008) , we establish 
a mass-dependent fit for M 2 oo c /M 2 oom 


M 2 00c , r , M 20 OC 

where 7 and <5 depend on $7 m and redshift as 


7(fl m , z) = 70 + 71 exp - 


72 - z 

73 


and 


(A2) 


<5(Om,X) 

= So + 81 Z 


(A3) 

70 (O m ) = 

3.54 x 1CT 2 

+ n° 09 


7l (O m ) — 

4.56 x 1CT 2 

+ 2.68 x io _2 /o m 


72 (^m) = 

0.721+ 3.5C 

1 x 10' 2 /O m 


73 (Om) — 

0.628+ O.164/O m 


5o(O m ) = 

—1.67 x 10~ 

" 2 + 2.18 x 10 _2 O m 


5i(n m ) = 

6.52 x 10" 3 

- 6.86 x 10“ 3 O m . 

(A4) 

accurate at 

the few percent level in the range 0 < 

z < 2, 


10 13 < M 200 c /M q < 2 x 10 16 , and 0.15 < ft m < 0.5. 


APPENDIX B: IMPACT OF INITIAL REDSHIFT 

The simulations used in this work are initialized at redshift zm = 
60 using the Zel’dovich approximation. Under this assumption, 
running a simulation from too low an initial redshift can lead to 
the suppression of the formation of high-mass haloes (|Reed et al.| 
| 20 T 3 ] >. 


We confirm that our simulations are converged by running 
two additional control simulations. These correspond to Boxl/mr 
DMonly as used in the main body of this work, except that we 
choose higher initial redshifts, z m j = 80 and z^ = 100. In Fig. m 
we compare the abundance of haloes in the three simulations at red¬ 
shifts z = 0.5 and z = 0. All three simulations agree very well 
within the error bars. We conclude that our simulations initialized 
at redshift Zi„i = 60 are converged and suitable for calibrating the 
HMF. 


APPENDIX C: SUBSTRUCTURE AT A 500 c 

As discussed in Section |4.2| the Tinker08 HMF predicts more 
haloes for SO A 200 C and Asooc at low masses M < 10 14 Mq than 
our DMonly fit, and both HMF agree much better for A 2 oom- This 
is mainly due to different assumptions about halo definition and 
identification. In this section, we show that our simulations are in 
good agreement with the low-mass end of the Tinker08 HMF at 
Asooc when haloes are extracted assuming the same halo definition 
as applied in Tinker08. The different definitions are: 

• In this work, we extract spherical halo masses using the min¬ 
imum in the gravitational potential within each halo as the center. 
In this approach, a halo with a given mass for one SO definition 
has exactly one counterpart in any other SO definition. Another 
consequence is that no subhalo that would actually exceed large 
overdensity thresholds like Asooc is extracted within a halo. 

• In Tinker08 SO masses are computed for all (potentially over¬ 
lapping) haloes whose centers are separated by more than r a • As 
a consequence, for high overdensities, and within massive haloes, 
particles in the main halo and in its subhaloes may be counted mul¬ 
tiple times. Furthermore, a halo identified at some SO can corre¬ 
spond to multiple exposed haloes at higher overdensity. 

We investigate the impact of substructure at Asooc in our sim¬ 
ulations by extracting another sample of haloes adopting the Tin- 
ker08 approach. We apply a SO finder to our DMonly simula¬ 
tions at redshift z = 0 and extract all haloes that exceed the SO 
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Asooc; this sample also includes subhaloes. In Fig. |B 1| we com¬ 
pare the abundance of objects extracted in this way to the abun¬ 
dance of main haloes as extracted by SUBFIND and used throughout 
this work. Both agree for masses above ~ 1O 13 M0, but at lower 
masses, subhaloes contribute to the abundance by up to ~ 25% 
for M 500 C — 10 12 Mq. This observation is in agreement with the 
discussion in Tinker08 and Kravtsov et al. (2004). Indeed, Fig. |Bl| 
suggests that, for masses below ~ 10 13 Mq, the total abundance of 
haloes including subhaloes extracted from our simulations agrees 
with the Tinker08 HMF. 

We conclude that, in the low-mass regime, the differences be¬ 
tween the HMFs presented in this work and the HMF by Tinker08 
are due to different assumptions on the halo definition. We note 
that these differences do not affect the shape of the HMF at higher 
masses M > 10 14 Mq which are important for the cosmologi¬ 
cal implications discussed in this work. However, this also makes 
it clear that the choice of (sub)halo identification applied to any 
observed cluster sample must be consistent with the identification 
method applied to the simulation data. For the cluster samples ob¬ 
served with the SPT and Planck, these differences are negligible, 
but they may indeed be important when using groups selected by, 
for example, eROSITA. 
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